clear variables; close all;

basename='Zoetwatervat';
save basename basename

fprintf('Zoetwater-opslagvat\n');
%% Get Variables
load setup
load CurrentSeriesNumber.mat;
n=CurNum;   % set current serie index


xSand =[    % [cm] Location sandlayer meausured with ginput on the pre scaled photo of the setup
   -0.0830
   97.5496
   97.1527
    0.0162
    ];
zSand =[
0.0806
-0.0186
16.7495
16.4519    
    ];



xWall =[    % Wall location for the Horizontal Flow Package
    SETUP_XWALL(n)
    SETUP_XWALL(n)
    ];

zWall =[
    SETUP_ZWALL(n)
    30      % [cm] 30 is top sandlayer, z-ax is measured with 0 [cm] the top of the impermeable layer
    ];   

%% The grid, the box is 96x96 cm and 1.8 cm thick
FRESHWATER=0;       % Relative minimum concentration
SEAWATER  =1;       % Relative maximum concentration
k=CONDUC(n)/(24*60);% CONDUC[cm/d]->[cm/min]
peff=POREFF(n)/100; % [] calibrated

MW   =98;           % [cm] Width of model
MH   =17;           % [cm]
D    = 10;          % [cm] thickness of model
Fringe=1;           % [m] Capillary fringe see page 34
zIface=17;          % [cm]

xzMilkSupply=[1,MH];% point of milksupply

%% Grid is always 3D, scale actual model to picture
dx=1.0;             % [cm] grid cell width
dy=  D;             % [cm] grid cell length = thickness of model; The model has become a 2D model, Hele Shaw
dz=1.0;             % [cm] grid cell size vertical

xGr=[0:dx:MW MW];	% [cm] grid, making sure that MW is always included irrespective of dx !
yGr=[0 D];          % [cm] grid
zGr=[MH:-dz:0 0];	% [cm] grid, making sure that 0 is always include irrespective of dz !

[xGr,yGr,zGr,xm,ym,zm,DX,DY,DZ,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);   % define gridproperties

[XM,ZM]=meshgrid(xm,zm);

[ixMilk,izMilk]=xyzindex(xzMilkSupply(1),xzMilkSupply(2),xGr,zGr);
[ixWall,izWall]=xyzindex(xWall(1),zWall(1),xGr,zGr);

Z=zeros(Ny,Nx,Nz);
for i=1:length(zGr)
    Z(:,:,i)=zGr(i);
end

%% Set initial concentration
if REUSECONC(n)==1, % if reuseconc==1 the final concentration of the previous cycle is beeing used
    load ConcSaved;
    STCONC=ConcentrationSaved;
    STCONC=permute(STCONC,[3,2,1]);
else
    STCONC=ones(Ny,Nx,Nz)*SEAWATER; % start concentration is all saltwater
end

%% Model arrays located on the picture, cell properties are being added through IBOUND matrix, see 
hCanL=MH; hCanR=15; % approximate canal water elevation (see workbook sheet PER)

IBOUND=ones(Nz,Nx);
% xLBnd=ixWall+1;   Not being used
% xRBnd=Nx;
IBOUND(Nz,Nx,1)=2;  % CHD package
% IBOUND(izMilk,ixMilk)=5;         % WELL package % milk injection point/Cells to which Well package applies
IBOUND(1,1:ixWall-1)=6;  % RECH package
% IBOUND(izFSSE,ixFSSE)=7;  % FSSE

%% Horizontal Flow Barrier
ILay=1; % Layer index
c=999999;
HFB=mf_setHFB([xWall,zWall],xGr,zGr,ILay,c);


%% Conductivity
HK = ones(size(IBOUND))*k; HK(IBOUND==0)=0;
VK=HK;
HK(HK>0 & ZM>hCanL+Fringe)=k/10; % HK=0 above full capillary zone

PEFF=ones(size(IBOUND))*peff;    % set peff for cells
PEFF(ZM>hCanR)=peff/3;           % unsaturated
PEFF(IBOUND==0)=0;             % inactive

% Permute Arrays to 3D
IBOUND=permute(IBOUND,[3,2,1]);
HK    =permute(HK    ,[3,2,1]);
VK    =permute(VK    ,[3,2,1]);
PEFF  =permute(PEFF  ,[3,2,1]);

ICBUND=ones(size(IBOUND));

STRTHD=ones(size(IBOUND))*-30;
STRTHD(Ny,Nx,Nz)=0;


[PerNams,PerVals]=getPeriods(basename);  % get stress periods from workbook
NPER = size(PerVals,1);                     % total number of stress periods
%% CHD boundaries for MODFLOW and PNTSRC for MTRDMS/SSM
[PerNams,PerVals]=getPeriods(basename);  % get stress periods from workbook
QMilk=[PerVals(:,strmatchi('QMilk',PerNams))];
% QMilk=[PerVals(:,strmatchi('QMilk',PerNams))]./PerVals(:,strmatchi('PERLEN',PerNams)); %QMilk defined in PER-Excel

IL=cellIndices(find(IBOUND==2),size(IBOUND),'LRC'); uL=ones(size(IL(:,1)));
IS=cellIndices(find(IBOUND==5),size(IBOUND),'LRC'); uS=ones(size(IS(:,1)));

ITYPE=2; % well
ITYPE2=-1; % CHD

iPer=1;

CHDOPT=1;      
CHDDENSOPT=1;  % Activate special CHD options in Seawat
CHDDEN=1029; % Density salt water[kg/m3]
CHDEN=1029; % Density salt water[kg/m3]
CHD=[ iPer*uL IL uL*[MH MH CHDDENSOPT CHDDEN]];

PNTSRC=[[iPer*uS IS uS*[FRESHWATER ITYPE FRESHWATER]];[iPer*uL IL uL*[SEAWATER ITYPE2 SEAWATER]]]; % milk injection point

%% WELL package
% [wel,WEL,PNTSRC,NPER]=mf_setwells(basename,xGr,yGr,zGr,HK,'wells');
% PNTSRC=[[iPer*uS IS uS*[FRESHWATER ITYPE FRESHWATER]];[iPer*uL IL uL*[SEAWATER ITYPE2 SEAWATER]]];

% NPER=size(PerVals,1);
% hL   =PerVals(:,strmatchi('hL',   PerNams));
% WEL=[ iPer*uS IS uS*QMilk(iPer) ]; % Milk injection point

% for iPer=2:NPER
%     if   hL(iPer) == hL(iPer-1), CHD=[CHD; -iPer NaN(1,5) CHDDENSOPT CHDDEN];
%     elseif hL(iPer)>0,           CHD=[CHD;  iPer*uL IL uL*[MH MH CHDDENSOPT CHDDEN]];
% 
%     end
%     if  QMilk(iPer)==QMilk(iPer-1), WEL=[WEL; -iPer NaN(1,4)];
%     elseif abs(QMilk(iPer))>0,      WEL=[WEL;  iPer*uS IS uS*QMilk(iPer)];
%     end
%     
%     PNTSRC=[PNTSRC; iPer*uS IS uS*[FRESHWATER ITYPE FRESHWATER]];
% end

% %% RCH Recharge see below N=4.2 L/h according to Pennink
% 
% W=45; % [cm] width of rain added to top of model
% M=37; % [cm] center of rain added to top of model (mid between 2 screws see photo)

%% RECH package [cm/IPER]->[cm/min]
NrMilk=length(find(IBOUND==6)); % amount of x cells where rech applies
N=QMilk*MULTIMILK(n)*(1/(dx*NrMilk))*(1/dy);  % recharge [cm/min] QMilk [cm3]

RECH=zeros(Ny,Nx,NPER);         % Allocate

for i=1:NPER
    RECH(:,IBOUND==6,i)=N(i);
end

%% save variables needed for mf_analyze

save underneath xSand zSand ixWall izWall dy dx dz ixMilk izMilk QMilk uS peff xGr yGr zGr